Delocalizing effect of the Hubbard repulsion for electrons 
on a two-dimensional disordered lattice 
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We study numerically the ground-state properties of the repulsive Hubbard model for spin-1/2 
electrons on two-dimensional lattices with disordered on-site energies. The projector quantum Monte 
Carlo method is used to obtain very accurate values of the ground-state charge density distributions 
with N p and N p + 1 particles. The difference in these charge densities allows us to study the 
localization properties of an added particle. The results obtained at quarter-filling on finite clusters 
show that the Hubbard repulsion has a strong delocalizing effect on the electrons in disordered 2D 
lattices. However, numerical restrictions do not allow us to reach a definite conclusion about the 
existence of a metal-insulator transition in the thermodynamic limit in two-dimensions. 

PACS numbers: 71.10.Fd, 73.20.Fz, 71.23.An 
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I. INTRODUCTION 

The interplay between disorder and electron-electron 
interactions has been a subject of intense activity in the 
last few years. Recent experiments have shown the ex- 
istence of an apparent metal-insulator transition (MIT) 
in two-dimensional (2D) semiconductor devices [1]. This 
observation came as a surprise to the community, since 
the scaling theory of localization, developed for disor- 
dered, non-interacting systems predicts insulating states 
even for infinitesimal disorder strengths [2]. Therefore, 
these experiments have promoted intense theoretical ac- 
tivity. There is at present no consensus and considerable 
controversy surrounds this problem [1]. A novel feature 
of the experimental high mobility samples investigated is 
their exceptionally low electronic density n s . This leads 
to an unusually high value of the dimensionless param- 

— 1/2 

eter r s oc n s of up to 80. The parameter r s sets the 
scale of electron-electron interaction E e _ e as compared to 
the Fermi energy Ep through r s k, E e _ e /Ep. Thus the 
question of electron-electron interaction effects becomes 
important in these disordered systems. 

The analytical treatment of the problem of disordered, 
interacting electrons is possible only in some limiting 
cases. The effects of weak interactions in disordered sys- 
tems have been studied in great detail in the metallic (dc- 
localizcd) regime [3]. The other extreme, corresponding 
to the strongly localized system can be treated by mean- 
field methods [4]. Although these analytical approaches 
have provided many useful physical results, the general 
treatment of disordered quantum many-body systems re- 
mains an unsolved problem. Indeed, one of the most 
powerful methods developed for the analytical treatment 
of disordered systems, namely, the supersymmetry ap- 
proach [5], cannot handle the effects of electron-electron 
interactions. Given this context, numerical approaches 
play a crucial role in the treatment of disordered, inter- 



acting systems. 

Several numerical approaches have been applied to 
the study of 2D disordered, strongly correlated systems 
[6-13]. Among the approximate approaches, Hartree- 
Fock calculations with residual interactions similar to the 
configuration interaction approaches of quantum chem- 
istry have been applied to systems with spinlcss fcrmions 
and electrons with spin [6,8]. Exact diagonalization ap- 
proaches have been useful but suffer from severe limita- 
tions in the accessible system size and number of particles 
[7,10,11,13]. From these studies, it has been observed 
that repulsive electron-electron interactions can have a 
delocalizing effect in small systems. In addition, it has 
been shown experimentally [1] that the spin degrees of 
freedom play a crucial role in the physics of these strongly 
interacting systems. The inclusion of the spin degrees 
of freedom renders the numerical calculations even more 
difficult and strongly reduces the number of fermions ac- 
cessible (see e.g. [6,10-13]). 

Quantum Monte Carlo (QMC) approaches provide a 
powerful alternative to the treatment of quantum many- 
body systems. These methods are in principle exact, 
apart from statistical errors and allow the treatment of 
much larger system sizes with many particles. The finite- 
temperature determinantal QMC approach has been ap- 
plied to the two-dimensional disordered Hubbard model 
and signatures of a metal-insulator transition were ob- 
tained [14]. However, these calculations were carried out 
at finite temperature and could not access the ground 
state of the system. In previous work, we studied the 
ground state properties of the disordered 2D Hubbard 
model by the projector quantum Monte Carlo (PQMC) 
method [15]. We studied the properties of the Green's 
function, charge density and inverse participation ratio 
against model parameters and system size. While we ob- 
served some local charge reorganization, we could not de- 
tect any significant delocalizing influence of the Hubbard 
repulsion U on the many-body ground state. However, it 
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should be noted that all the physical quantities studied 
in Ref. [15] were obtained by integrating over all parti- 
cles and effectively, the entire energy spectrum. Here we 
introduce a new approach which allows us to study the 
properties of a single added particle and therefore empha- 
sizes the physical effects of interactions in the proximity 
of the Fermi edge. 

This approach uses the inverse participation ratio 
(IPR) extracted from the charge density differences of 
two many-body ground states, as described below. The 
IPR, £, for a normalized single-particle wavefunction ip(i) 
is given by £ = [^2 i where i is the site index of 

the system. Clearly, VKOI 2 can be identified as the one- 
particle charge density at the site i. This definition is usu- 
ally carried over to many-body systems by renormalizing 
the total charge density at site i, p(i,N p ), which is ob- 
tained in the standard way from the ground state many- 
body wavefunction for N p particles. The IPR £ is then 
calculated with the effective charge density p(i,N p )/N p . 
We found that the IPR obtained from such a definition 
showed slight variations with model parameters and prac- 
tically no variation with system size [15]. We believe 
that the physical reason for this weak variance is the fact 
that this procedure effectively integrates over all ener- 
gies and therefore the dominant contribution comes from 
the states deeply below the Fermi energy, which remain 
strongly localized even in the presence of interactions. 
Therefore, it is necessary to find a method which is more 
sensitive to the contribution of states in the vicinity of 
the Fermi energy. Recently, we studied the localization 
properties of the disordered, attractive Hubbard model in 
two- and three-dimensions [16]. In Ref. [16], we showed 
that an IPR calculated for an added pair of particles is a 
very sensitive and relevant quantity to study the localiza- 
tion properties of the wavefunction. Based on our results 
for the attractive Hubbard model, we have introduced 
a related quantity in the repulsive case, the IPR for a 
single added particle, which turns out to be much more 
sensitive to the effects of interaction. This quantity has 
certain similarities to the single-particle tunneling am- 
plitude. The latter has recently been studied by exact 
diagonalization methods for small clusters with spin- 1/2 
fermions [13]. This work provided indications for a delo- 
calizing effect induced by repulsive interactions in disor- 
dered systems. However, the number of particles studied 
was rather restricted due to the limitations of the exact 
diagonalization methods. With our method and exten- 
sive, highly accurate PQMC simulations, we study con- 
siderably larger numbers of particles in presence of strong 
electron-electron interactions. In this study we report a 
significant delocalizing influence of the Hubbard repul- 
sion on 2D disordered electronic systems. This paper is 
organized as follows : after this Introduction, we describe 
the method used and the tests performed in the next sec- 
tion. Our results and discussion are presented in detail 
in the third section. 



II. MODEL AND METHOD 

The Hamiltonian studied in this paper is the disor- 
dered, repulsive Hubbard model given by: 

H = H A + Hj 
where c\ a (ci a ) creates (destroys) an electron at site i 



with spin a and n. 



is the corresponding oc- 



cupation number operator. The hopping term t between 
nearest neighbor lattice sites characterizes the kinetic en- 
ergy and the random site energies ti are taken from a box 
distribution over [— W/2,W/2]. The parameter U mea- 
sures the strength of the screened, repulsive Hubbard in- 
teraction (U > 0). We have considered both the one- 
and two-dimensional cases, with periodic boundary con- 
ditions in all directions. In the 2D case, the sites i lie 
on a rectangular lattice of linear dimension N x ,N y . The 
system size ./V = N x x N y then follows accordingly in one- 
and two-dimensions. In the limit U — 0, this Hamilto- 
nian reduces to the Anderson model (given by Ha) which 
is a standard model for the study of disordered systems 
[17]. In the absence of disorder, W = 0, this Hamilto- 
nian reduces to the usual Hubbard model, which is one of 
the best-studied models for correlated electronic systems 
[18]- 

We have studied this Hamiltonian by the projector 
quantum Monte Carlo (PQMC) method [19] and exact 
diagonalization calculations. The PQMC method was 
initially developed to study the ground state of the Hub- 
bard model (clean limit of Eqn. (1)). The method can be 
generalized, in principle quite simply, to incorporate dis- 
order via the random site energies. However, the actual 
implementation and convergence of the algorithm [15] is 
highly non-trivial compared to the pure case, as will be 
discussed in detail below. 

The PQMC method consists in filtering out the true 
ground state |^o) of the many-body system from an ap- 
propriately chosen trial function \<f>): 



1 ' e^cx> ^/(4>\e- 2eH \ 



(2) 



This method is exact in principle, apart from statistical 
errors and the sign problem which appears for fermions 
at U > 0. The Hamiltonian plays the role of the projec- 
tion operator through the term e~ eH , where plays the 
role of the projection parameter. The trial wave-function 
is usually formed as a product of up and down spin 
states from the eigenstates of the non-interacting Hamil- 
tonian. In our case, we chose the Fermi sea of Ha as the 
trial wave-function. In the PQMC procedure, the projec- 
tion operator exp(— OH) is first Trotter decomposed as 

(exp(-ArH A ) exp(-Ar.Hi)) with 6 = At x L. This 
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introduces a systematic error of order (At) 2 due to non- 
commutation of Ha and Hj . We have used the symmet- 
ric Trotter decomposition, which introduces a systematic 
error of (At) 3 . The interaction is then decoupled by 
a discrete Hubbard-Stratonovich transformation, by the 
introduction of TV x L Ising-like fields. This Ising model 
with complicated effective interactions is then treated by 
a Monte Carlo (MC) procedure to obtain the ground 
state properties of the system. The quantity calculated 
during the simulation is the zero-temperature, equal-time 
Green's function, dj = ^2 a (ipo\cl a Cj^\ipo) ■ During the 
simulation this Green's function can be used to obtain 
all the other static correlation functions. In the algo- 
rithm used 0(N 2 ) operations are required to update the 
Green's function after a MC step. This procedure intro- 
duces cumulative errors and therefore the Green's func- 
tion has to be recalculated from scratch regularly during 
the simulation (every Lc steps), which requires 0(N 3 ) 
operations. When the projection parameter becomes 
large, which is necessary for good convergence, the dif- 
ferent components of the wavefunction tend to become 
parallel during the projection process. Therefore, it is 
necessary to reorthogonalise the components of the wave- 
function regularly, every Lr time steps. In addition, it 
is well known that quantum simulations of fermionic sys- 
tems suffer from the sign problem except in some special 
cases. However, it has been observed that disorder in 
fact diminishes the magnitude of the sign problem. In 
our simulations, the sign problem is well under control, 
with the number of negative signs being less than 1% of 
the total number of steps considered. 

We have studied systems with particle number (N p ) up 
to 25 fermions on lattice sizes of up to 6 x 8 sites, with 
particle density always around quarter-filling. The sim- 
ulations were carried out in the S z = sector for even 
numbers of particles and the S z = 1/2 sector for odd 
number of particles. The specific quantity studied in this 
paper is the difference of charge densities Sp(i) between 
systems with N p and N p + 1 particles (see the next sec- 
tion for more details). Here, we discuss the convergence 
of the PQMC calculations. The charge densities involved 
in the calculation 5p(i) are obtained from two indepen- 
dent simulations of the same disorder realization with 
different particle numbers, N p and N p + 1. Therefore, it 
becomes necessary to measure very accurately the distri- 
bution of this added particle over the lattice. Clearly, it 
is harder to measure 5p(i) accurately than measuring the 
total charge densities p(i, N p ) in the ground state with 
N p particles. We carried out extensive tests to verify the 
quality of our 8p(i) data, by varying the PQMC param- 
eters until convergence was obtained. We have tested for 
convergence by varying the following PQMC parameters: 
the projection parameter 0, the Trotter time-step At, 
the Monte Carlo parameters, the reorthogonalization in- 
terval Lr and the interval to recalculate the Green's func- 
tion Lc- We have gone up to = 15 and At = 0.05. For 
the physical parameters used in this paper, we find that 



= 10 with At = 0.1 are sufficient. This corresponds to 
a systematic error of 10~ 3 . After varying the parameters 
Lr and Lc we decided to recalculate the Green's function 
with reorthogonalization of the components, after every 
5 time-steps ( Lr = Lc = 5). We have also checked the 
MC parameters and find that 3000 sweeps are adequate 
for convergence, with 1000 sweeps for equilibration. 

The disorder average was carried out over Nr different 
disorder realizations with Nr = 16 in the PQMC simu- 
lations and Nr = 100 for most exact calculations. The 
site-energies are randomly chosen for Nr disorder real- 
izations at W/t = 1 and then scaled proportionally to 
W/t for stronger values of disorder. 

While our main motivation is the study of the 2D case, 
there are many useful reasons to study ID systems. The 
ID case can be studied conveniently by exact diagonal- 
ization methods, while changing the system size. Thus, 
we can study the variation of various properties such as 
the inverse participation ratio (IPR) £, with system size, 
exactly. These calculations provide a strong independent 
check on the PQMC data at small sizes. Further, the lo- 
calization effect is much stronger in ID and Sp is usually 
strongly peaked in ID as compared to 2D. Therefore, if 
the algorithm is capable of reproducing a localized peak, 
this provides a strong check on the method, at the given 
range of physical parameters (examples provided in the 
next section). 

III. RESULTS AND DISCUSSION 

In order to study the localization properties of the 
system, we use the charge density distribution for an 
added particle at the Fermi level, given by: Sp(i) = 
p{i,N p + 1) — p(i,N p ), where p(i) is the ground state 
charge density at site i. The values of p(i, N p ), p(i, N p +1) 
are obtained from two independent PQMC simulations 
for the same disorder realization. At U = 0, this 8p(i) is 
identically equal to the one-particle probability distribu- 
tion (added particle) at the Fermi edge. In the interacting 
case, this charge density distribution (8p(i)) is not gen- 
erally equal to a probability distribution. For example, 
it is not necessarily positive definite. However, in all the 
cases studied we find that this quantity is always positive. 
Therefore, 8p(i) can be considered to be similar to a one- 
particle probability distribution. Thus, we can associate 
an inverse participation ratio (IPR) for an added parti- 
cle, for a given disorder realization as: £ = (yv 5p(i) ) . 
The disorder averaged IPR is given by (£). At U = 0, 
the IPR is a standard tool used to obtain the number 
of sites over which the particle is localized [5] . We have 
successfully used this approach in previous work on the 
disordered, attractive Hubbard model for a quantitative 
description of the localization properties of the ground 
state [16]. 
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FIG. 1. IPR (£) versus interaction strength U/t, at 
W/t — 7, for N p = 6 particles on a one-dimensional lat- 
tice with N = 12 sites (triangles, average over 100 disorder 
realizations), and on a 4x3 lattice (circles, average over 16 dis- 
order realizations). Data comes from exact diagonalization. 
Here and in the following figures error bars denote statistical 
errors. 




FIG. 2. IPR {£) versus disorder strength W/t. Continuous 
curves are for the ID system of 12 sites and N p — 6 particles, 
at U/t = (open circles) and U/t = 4 (filled squares) . Curves 
with dashed lines are for the 4x3 lattice and N p = 6 particles, 
at U/t — (open circles) and U/t — 4 (filled squares). Data 
come from exact diagonalization and are averaged over 100 
disorder realizations (except the point for the 4x3 lattice at 
U/t = 4, W/t = 12, with N R = 16). The continuous curve 
with inverted triangles is the ratio {£(U/t))/{£(U/t = 0)), at 
U/t = 4 for ID and the dashed curve with triangles the same 
ratio for 2D. 

In Figs. 1-2, we study the behavior of the IPR as a 
function of U/t and W/t to try and establish the inter- 
esting range of physical parameters in this system. In 



Fig. 1, we compare the IPR (£) for a ID ring of 12 
sites and a 4 x 3 lattice in 2D as a function of interac- 
tion strength U/t. From the figure we see that increasing 
U /t from tends to increase the IPR up to intermediate 
strengths of the interaction. Thus, the optimal value ap- 
pears to be around U/t = 4 in both ID and 2D. Indeed, 
for U/t > 4, the value of IPR starts to decrease. In Fig. 
2, we see the variation of (£) for small system sizes for 
5 < W/t < 12. Since the ratio (£(U/t = 4)>/(£(t//i = 0)) 
is practically constant, this justifies our choice of param- 
eter range for W/t. Disorder strength W/t < 5 would 
lead to states with localization length larger than the ac- 
cessible system size even at U/t = 0. Convergence of the 
PQMC calculations becomes progressively more difficult 
for larger disorder strengths. Therefore, we studied the 
variation of the IPR for W/t = 5, 7 in the 2D case. 




FIG. 3. Charge density difference 6p(i) for a ID lattice with 
N = 12, N p = 6, W/t = 7, U/t = (thin line, inverse partic- 
ipation ratio £ = 4.31), exact result for U/t — 2 (thick solid 
line, £ = 6.58), and quantum Monte Carlo data for U/t = 2 
(dashed line, £ = 6.93). 

Thus, it would seem ideal to study the IPR of ID and 
2D systems for 5 < W/t < 7 and U/t « 4. In ID we 
could access these parameter values conveniently for the 
system sizes studied by exact calculations. However, this 
becomes much more difficult in 2D. This is because we 
study a small difference of large total charge densities in 
the PQMC simulations. We found that the accuracy of 
our method for Sp(i) is good for for values of the inter- 
action U/t < 2. While we are not exactly at the optimal 
value of U/t, we are nevertheless in a region with suffi- 
ciently strong interactions. Indeed, it can be seen from 
Fig. 1 that U/t — 2 already has a substantial delocalizing 
influence on the system. In Fig. 3, we present the PQMC 
calculation of Sp(i) in ID at U/t = 2 as compared to ex- 
act calculations. We note that the Sp(U = 0) is strongly 
peaked. Increasing U/t to 2 radically changes the picture 
and shifts the peak completely. The PQMC curve shown 
reproduces the quantitative picture of charge density dif- 
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ference. It should be noted that the calculation begins 
with a trial wavefunction corresponding to the U/t = 
data and changes completely to give the correct physical 
picture. Quantitatively, the usual errors seen at these 
values of the physical parameters were around 3 — 5%. 
The 2D case is more delocalized compared to ID and we 
have observed that convergence is better in 2D. For ex- 
ample, we have (£) = 7.19 (exact) and 7.43 (QMC) for 
a 4 x 3 lattice at U/t — 2 and W/t = 7, averaged over 
16 disorder realizations. This corresponds to an error of 
about 3% for the most extreme parameter values studied. 
However, for U/t = 4, the error increases to 8 — 10% and 
therefore we restrict our studies to U/t < 2 in 2D. 
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FIG. 4. Charge density difference Sp for a ID lattice with 
N = 12, N p = 6, W/t = 7, and U/t = 4: p{N p + 1) - p{N p ) 
(continuous curve with circles) and p(N p + 2) — p(N p + 1) 
(dot-dashed curve with squares). Data are from exact diago- 
nalizations, for the same disorder realization as in Fig. 3. For 
comparison, see data for Sp at U/t = in Fig. 3 (thin line). 

Since we study a model of electrons with spin, it is 
interesting to consider even-odd effects in the particle 
number N p and system size N. In Fig. 4, we consider 
the effect of progressively adding one and two particles 
to a 12 site ring with 6 electrons. From Fig. 3 and Fig. 
4, we see that the first added electron at U/t = 2 and 
4, has an entirely different peak as compared to the non- 
interacting case. In Fig. 4, we see that the second added 
electron also occupies a different region, as seen from the 
position of the peak. This leads to the following physical 
image for the ID case : each added electron finds its 
own location in space with optimal energy and it is well 
localized by external disorder and random distribution of 
charges of other electrons. This is in contrast with the 
non-interacting case, where the non-interacting orbital 
remains the same for odd and next even electron. 

The situation is qualitatively different in 2D. Indeed, in 
Fig. 5, we show the charge density difference for one and 
two added particles on a 6 x 6 lattice with 18 electrons. 
This quantity is obtained from PQMC simulations of the 



system for a particular disorder realization, taken at W/t 
= 7. The data show that the initial U/t — configuration 
is very clearly localized for the given disorder realization 
and lattice size. It is seen from the figure that the in- 
troduction of a repulsive Hubbard interaction (U/t = 2) 
leads to a substantial derealization of the added parti- 
cle. This is also borne out quantitatively, since the IPR 
increases practically by a factor of 3, as compared to 
U/t = 0. In the non-interacting case, both added parti- 
cles occupy the same orbital. Therefore, the peak for the 
second added particle, is trivially identical to the first. 
It is remarkable that in the interacting case, the peak is 
transformed to a much more extended distribution over 
the lattice. Furthermore, the second added particle has 
practically the same distribution Sp(i) as the first. Thus, 
there appears to be almost zero effective repulsion be- 
tween the two added particles, despite the interaction 
strength U/t = 2. This is completely in contrast to the 
scenario in ID with interactions, where we observe sig- 
nificant repulsion between the two added particles. We 
note that this is the case for every disorder realization 
studied. 




FIG. 5. Charge density difference 8p for a 6 x 6 lattice 
with N p = 18 particles, W/t = 7, U/t = (left) and U/t = 2 
(right): p{N p + 1) - p(N p ) (top) and p(N p + 2) - p(N p + 1) 
(bottom). Inverse participation ratios are £ = 6.2 (top left), 
I = 18.3 (top right), £ = 6.2 (bottom left), and £ = 16.7 
(bottom right). 

We now turn to a more quantitative picture of the one 
and two dimensional systems. For this, it is necessary 
to average the data over different disorder realizations 
and to vary the system size. We start the analysis from 
the data for the ID case. In Fig. 6, we present the 
IPR for interacting particles on ID rings of 4-12 sites at 
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constant density of particles (quarter filling) and compare 
against the variation of the IPR at U/t = 0. We find a 
derealization effect even for ID rings. We note that the 
curve for U/t = is well saturated in ID even at 8- 
12 sites. The U/t = 2 curve shows signs of saturation. 
The derealization effect is most visible for U/t = 4, as 
expected from Fig. 1, where unfortunately, we do not 
have access to PQMC data for comparisons. The ratio 
(Z(U/t))/(£(U/t = 0)) goes to 1.64 for U/t = 2 and 1.94 
for U/t = 4. 




FIG. 6. Inverse participation ratio (£) versus system size 
N, for ID chains with periodic boundary conditions, at 
W/t = 7, U/t = (circles), U/t = 2 (diamonds), and U/t = 4 
(squares). Triangles and inverted triangles give the ratio 
k(U/t))/(£(U/t = 0)), at U/t = 2 and U/t = 4, respectively. 
Data come from exact diagonalization and are averaged over 
100 disorder realizations. 

In 2D previous exact diagonalization studies have gone 
up to 6 electrons on a 4 x 4 lattice, or 4 electrons on a 
6x6 lattice [10,11,13]. Therefore, it is of great impor- 
tance to access larger system sizes with more particles. 
In Fig. 7, we present the IPR for different lattice sizes 
in 2D, obtained from PQMC simulations. We have gone 
up to 48 sites (6x8 lattice) and 24-25 particles, which 
goes beyond any existing study of the ground state in 
the literature. We observe a dclocalizing effect of the 
Hubbard repulsion manifested by an increase of the IPR 
with system size and interaction strength. There is a re- 
markable difference between the present result and the 
IPR for the full many-body ground state wavefunction 
obtained in Ref. [15]. In the previous work, we could 
notice no size effects over a large range of W and lattice 
sizes. In fact, the curves for different lattice sizes against 
U/t at a given value of W/t tended to collapse as seen 
in Figs. 8 and 9 of Ref. [15]. On the contrary, in the 
present work, we see from Fig. 7 that effects of lattice 
size on the IPR for an added particle are considerable. 
Furthermore, we see that the Hubbard repulsion has a 



strong dclocalizing influence, compared to the data for 
the non-interacting case. As explained in the previous 
sections, we attribute the difference between the present 
results and those presented in Ref. [15] to the fact that 
the IPR from Sp measures directly the response in the 
vicinity of the Fermi level. 
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FIG. 7. Inverse participation ratio (£) as a function of the 
number of sites N, for a 2D lattice, with quarter filling (i.e., 
N p = 6, 12, 18, and 24 fermions on 4 x 3, 6 x 4, 6 x 6, and 
8x6 lattices, respectively), with U/t — (empty symbols, 
average over 10 3 disorder realizations), U/t = 2 (filled sym- 
bols, average over 16 disorder realizations), and W/t = 5 
(circles) and W/t = 7 (triangles). The inset shows the ratio 
k(U/t))/{£(U/t = 0)), at U/t = 2, for W/t = 5 (circles) and 
W/t = 7 (triangles). PQMC parameters: 6 = 10, At = 0.1, 
except for 8 x 6 lattices where O = 12, At — 0.08. 

The ratio (£(U/t))/(£(U/t = 0)) can be considered as 
a quantitative measure of derealization effect induced 
by repulsive interaction. From the inset of Fig. 7, it 
can be seen that this ratio rises with system size for a 
given density and we do not observe saturation at the 
system sizes studied. Our data for this ratio for small 
clusters in 2D are comparable to the enhancement ratio 
obtained in Ref. [13] though the quantity studied in Ref. 
[13] was slightly different i.e. the one-particle tunneling 
amplitude. From Figs. 6 and 7 we can compare the val- 
ues of the ratio (£(U/t))/(£{U/t = 0)) in ID and 2D. 
At the system sizes studied these values are comparable 
for U/t = 2. However, it would necessary to analyze the 
behavior of this ratio for larger system sizes in order to 
get a clear answer to whether the localization properties 
are indeed different in ID and 2D. At the same time, we 
note that there seem to be qualitative differences based 
on even-odd effects for added particles in ID and 2D (see 
Fig. 4 and Fig. 5 and discussion there). This differ- 
ence favors the picture of stronger derealization in 2D 
as compared to ID. 

Even if our data clearly show a repulsion induced der- 
ealization effect, they do not permit us to draw a definite 
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answer about the existence of a metal-insulator transi- 
tion for this system in the thermodynamic limit. Indeed, 
even though we have a significant number of fermions, we 
cannot go to larger system sizes because of the accuracy 
constraint of calculating charge differences. It should also 
be pointed out that our calculations done at U/t < 2 
are below the optimal value of the interaction strength. 
We expect that the delocalizing effect is even stronger at 
U/t = 4. 

In conclusion, we have studied the ground state prop- 
erties of the repulsive Hubbard model with disorder 
through the powerful PQMC method. Highly accurate 
simulations permit us to obtain the difference of charge 
density between two ground states with N p and N p + 1 
fermions respectively. The analysis of this characteristic 
clearly shows that the Hubbard repulsion has a delocaliz- 
ing effect in the system. We have observed some qualita- 
tive differences between ID and 2D for this characteristic. 
However, the restrictions on system size and interaction 
strength do not permit us to draw a definite conclusion 
about the existence of a MIT in the thermodynamic limit 
in 2D. 
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